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Abstract 

Solder  drops  spreading  on  metallic  substrates  are  a reactive  form  of  the  wetting  problem.  A 
metallic  component  may  diffuse  in  the  liquid  toward  a metal  substrate  where  it  is  consumed  by  a 
reaction  that  forms  a solid  intermetallic  phase.  The  resulting  spatial  variation  in  the  composition 
of  the  drop  may  cause  composition  gradients  along  the  free  surface  of  the  drop.  Together  with 
any  thermal  gradients  along  the  free  surface,  Marangoni  effects  may  in  turn  modify  the  bulk 
transport  in  the  spreading  drop.  Motivated  by  this  situation,  we  extend  lubrication  theory 
for  the  spreading  of  thin  drops  in  the  presence  of  gravity  and  thermocapillarity  to  include 
mass  transport  and  solutocapillarity.  We  use  an  approximate  solute  profile  in  the  drop  to 
derive  coupled  evolution  equations  for  the  free  surface  shape  and  concentration  field.  Numerical 
solutions  for  the  non-reactive  (single  component)  drop  agree  well  with  previous  theory.  In  the 
reactive  case,  we  are  only  able  to  compute  results  for  parameters  outside  of  the  range  for  solder 
materials.  Including  reactive  effects  in  the  model  impacts  the  flow  patterns  and  spreading  rates 
at  relatively  early  times;  but  by  the  end  of  the  spreading,  solutal  effects  have  died  out  in  the 
model. 


’Corresponding  author,  Tel.:  (301)  975-5913,  Fax:  (301)  990-4127 
^Technology  Administration,  U.S.  Department  of  Commerce,  Washington,  D.C. 
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1.  Introduction 


The  spreading  of  a droplet  on  a horizontal  surface  has  been  the  subject  of  much  research;  it  is  a 
relatively  simple  configuration  to  evaluate  the  wetting  of  a surface  by  a liquid.  Such  a configuration 
is  also  used  in  the  evaluation  of  wetting  in  solder  (or  braze)  systems,  where  a disc  of  solid  solder 
is  melted  on  a hot  substrate,  and  the  subsequent  rate  and  extent  of  the  spreading  is  measured.1-5 
While  there  is  not  a direct  correlation  between  the  wettability  of  the  solder  and  the  quality  of  the 
joint,  good  wetting  is  requisite  in  the  formation  of  a solder  joint  of  the  required  geometry.  In  solder 
systems  (e.g.  Pb-Sn  or  Bi-Sn  on  Cu  or  Au),  the  wetting  problem  is  complicated  by  a number  of 
effects,  including  heat  and  mass  transport  in  the  spreading  drop,  the  use  of  flux  to  remove  non- 
wetting oxide  patches  on  the  substrate,  as  well  as  reactions  in  which  intermet allies  form  at  the 
substrate  (e.g.,  Cu+Sn— >Cu6Sns).  The  liquid  surface  tension  on  the  free  surface  (as  well  as  other 
properties)  depend  strongly  on  the  temperature  and  composition  of  the  molten  solder,6-8  and  these 
variations  may  in  turn  affect  the  fluid  dynamics  of  the  spreading  through  the  generation  of  shear 
stresses  along  the  free  surface,  i.e. , the  Marangoni  effect. 

In  some  recent  experiments,  it  was  observed  that  the  contact  angle  in  the  final  shape  of  the 
solder  drop  was  in  the  range  of  2 to  16  degrees  (see  e.g.,  Ref.  1,  9).  At  these  small  angles,  the 
drop’s  height  is  much  smaller  than  its  radius,  and  so  there  is  a separation  of  scales;  the  lubrication 
approximation  may  be  applied  to  the  equations  of  fluid  motion  to  describe  the  spreading  of  the  thin 
drop.  The  presence  of  the  contact  line  for  spreading  drops  complicates  the  problem.  The  physics 
of  the  contact  line,  even  without  reactions  and  other  processes  to  complicate  matters,  is  an  area  of 
current  research.  In  mathematical  models  of  fluid  motion  that  use  the  no  slip  boundary  condition, 
there  is  a singularity  in  the  shear  stress  in  the  fluid  at  the  contact  line;10,11  this  singularity  may  be 
relieved  by  allowing  some  slip  of  the  fluid  along  the  substrate.12,13  Fortunately  for  the  macroscopic 
modeling  of  the  fluid  flow,  only  the  flow  very  near  the  contact  line  is  strongly  affected  by  the 
form  of  the  slip  applied  at  the  substrate12, 14  while  the  flow  far  from  the  contact  line  is  practically 
indistinguishable  for  different  slip  models. 

Though  the  contact  angles  for  the  spreading  drops  in  solder  experiments1, 2 are  small,  it  appears 
that  the  drops  only  partially  wet  the  substrate.  When  equilibrium  contact  angles  are  zero  (per- 
fect wetting),  a precursor  foot  is  predicted15  and  found  experimentally16,17  for  organics  on  inert 
substrates.  A precursor  foot  of  this  type,  with  a thickness  of  hundreds  of  Angstroms  or  less,  is 
not  observed  or  predicted  for  nonzero  values  of  the  equilibrium  contact  angle  (partial  wetting).  In 
the  perfect  wetting  case,  the  contact  line  position  is  given  by  a temporal  power  law;  the  long  time 
behavior  does  not  obey  power  law  spreading  in  the  partial  wetting  case. 

We  believe  that  the  Marangoni  effect  may  play  a role  in  the  spreading  of  solder  drops  because 
the  transport  of  heat  through  the  drop  sets  up  tangential  temperature  gradients  along  the  surface  of 
the  drop,  which  in  turn  cause  surface  tension  gradients  there.  These  gradients  cause  shear  stresses 
that  either  promote  or  retard  spreading.  Often  the  solder  drop  melts  and  spreads  on  a heated 
substrate,  tending  to  retard  spreading  for  most  solder  materials.  In  other  situations,  however,  the 
solder  spreads  on  a relatively  cool  substrate  while  heat  is  supplied  from  the  surroundings,  tending 
to  promote  spreading. 

Tangential  gradients  in  the  alloying  element  concentration  may  also  cause  surface  tension  gra- 
dients. Near  the  eutectic  point  for  Pb-Sn  solder  drops,  the  change  in  surface  tension  for  a 20K 
temperature  change  is  roughly  equal  to  a lwt%  concentration  change;6,7  this  suggests  that  the  con- 
centration variation  may  be  just  as  significant  as  temperature  change  with  regard  to  the  Marangoni 
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effect  in  the  reactive  spreading  of  solder  drops.  Implicit  in  this  argument  is  the  neglect  of  the  solder 
flux  which  is  often  used  in  soldering;  as  a first  model,  we  are  ignoring  the  role  of  the  solder  flux  in 
the  fluid  dynamics  of  the  spreading.  Our  thought  here  is  to  consider  Marangoni  effects  and  to  eval- 
uate their  relative  importance  during  the  spreading  process.  It  turns  out  that  the  real  parameters 
for  solder  systems  are  beyond  the  reach  of  the  theory  and  computational  method  developed  here, 
but  we  use  the  solder  spreading  process  as  motivation  for  considering  a model  of  reactive  spreading. 
As  much  as  possible,  we  base  the  parameters  we  use  on  the  Pb-Sn  solder  system. 

The  use  of  lubrication  theory  for  thin  drops  was  first  carried  out  by  Greenspan,18  and  has 
been  employed  by  many  workers  since.  In  Greenspan’s  work,18,19  an  evolution  equation  for  the 
shape  of  the  free  surface  of  the  drop  is  derived;  the  speed  of  the  contact  line  is  proportional  to 
the  difference  between  the  dynamic  contact  angle  and  the  advancing  equilibrium  contact  angle. 
In  the  limit  of  small  capillary  number,  that  evolution  equation  is  subsequently  solved  to  yield  an 
ordinary  differential  equation  for  the  contact  line  position.  In  a series  of  papers  by  Hocking  and 
coworkers, 20-22  the  contact  angle  is  assumed  constant  and  matched  asymptotic  expansions  in  the 
limit  of  small  slip  coefficient  are  used  to  derive  an  ordinary  differential  equation  for  the  contact 
line  position.  In  this  approach,  the  spreading  of  the  drop  is  limited  by  viscous  dissipation  and  the 
slip  at  the  substrate.  One  discussion  of  the  relationship  between  these  approaches  may  be  found 
in  Ref.  23.  Other  applications  of  lubrication  theory  in  spreading  include  Tanner,24  de  Gennes  and 
coworkers,15’25  and  Brochard-Wyart  et  al.26 

We  base  our  model  most  closely  on  that  developed  by  Ehrhard  and  Davis;27  they  suppose  that 
the  contact  line  mobility  limits  the  motion  of  the  drop  and  they  use  a modified  “Tanner  law”24 
as  a constitutive  relation  to  determine  the  motion  of  the  contact  line.  They  extend  previous  work 
to  include  a Marangoni  effect  due  to  temperature  gradients  in  the  drop.  Their  theory  has  been 
convincingly  verified  for  silicone  oil  on  glass  in  recent  experiments  carried  out  by  Ehrhard.28  The 
theory  has  been  recently  extended  by  Smith29  to  study  drop  migration  in  a horizontal  temperature 
gradient  and  by  Anderson  and  Davis30  for  volatile  droplets  on  a heated  plate.  Thermocapillary 
migration  has  also  been  studied  by  Ford  and  Nadim.31 

In  this  paper,  we  attempt  to  study  reactive  spreading  of  a drop  by  modeling  the  heat  and  mass 
transport  during  the  spreading  process,  and  include  a simple  model  for  the  reaction  at  the  substrate. 
We  generalize  the  work  of  Ehrhard  and  Davis27  to  include  gradients  in  the  concentration  (of  Sn,  say, 
in  Pb-Sn  drops)  caused  by  the  reaction  with  the  substrate  (Cu  and  Cu-Sn  alloys  in  Ref.  1,2).  We 
include  this  variation  by  deriving  a second  evolution  equation  for  the  alloying  element  distribution. 
In  order  to  derive  the  second  equation,  the  spatial  distribution  of  the  alloying  element  is  assumed 
to  be  a polynomial  in  the  vertical  coordinate;  that  is,  a Karman-Polhausen  type  approach  is  used 
(see  e.g.,  Ref.  32).  This  approach  was  motivated  by  the  study  of  spin  coating  with  an  evaporating 
solvent  carried  out  by  Reisfeld  et  al .33>34  The  resulting  coupled  evolution  equations  for  the  drop 
shape  and  alloy  distribution  are  solved  numerically  using  an  implementation  of  the  method  of 
lines;  our  numerical  solution  method  is  similar  to  that  used  by  Haley  and  Miksis35  for  a series  of 
drop  spreading  problems.  From  these  results,  and  from  work  in  progress,36  we  wish  to  assess  the 
importance  of  the  solutal  Marangoni  effect  in  reactive  spreading  of  drops. 

In  the  next  section  we  formulate  the  problem;  in  Section  3,  we  derive  a lubrication  theory  for 
the  reactive  spreading  of  a planar  (two-dimensional)  drop  on  a flat  substrate.  In  Section  4,  the 
numerical  methods  used  to  integrate  the  evolution  equations  are  discussed.  Section  5 gives  results 
for  reactive  and  non-reactive  spreading  drops;  discussion  and  conclusions  are  given  in  the  final 
section. 
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2.  Formulation 


We  consider  a two-dimensional  drop  in  a Cartesian  coordinate  system  ( x',z ')  on  a smooth,  hor- 
izontal plate  in  the  plane  z'  — 0.  The  primes  denote  dimensional  variables.  The  free  surface  of 
the  drop  is  given  by  z'  — h'(x',t')  and  the  contact  line  at  = 0 is  located  at  x'  — a'[t'). 

The  drop  is  assumed  to  be  symmetric  about  the  line  x'  — 0.  The  equations  governing  the  velocity, 
temperature,  and  concentration  fields  in  the  drop  are 


V'  V 

= 0, 

(2.1) 

— —Vp'  + /iV,2v'  - pgk, 

(2.2) 

= A V2T', 

(2.3) 

vv 

= dvY, 

(2.4) 

pcP 


where  k = (0, 1),  v'  = w ')  is  the  velocity  vector,  p'  is  the  pressure,  T'  is  the  temperature,  and  c' 

is  the  concentration  of  Sn  in  the  liquid.  Here  g is  the  magnitude  of  the  gravitational  acceleration, 
p is  the  density,  p is  the  dynamic  viscosity,  A is  the  thermal  conductivity  of  the  drop,  cp  is  the 
specific  heat,  and  D is  the  solute  diffusivity;  all  but  g are  properties  of  the  liquid. 

On  z'  — 0 (the  substrate),  we  have 


w - 0,  (5\h')u'zi  = v! , and  T'  — Tw. 


(2.5) 


The  slip  coefficient  /3'(h')  may  depend  on  b!  and  will  be  specified  later.  Our  Dirichlet  condition  on 
the  temperature  field  implies  that  the  substrate  is  a perfect  conductor. 

For  the  concentration  boundary  condition,  we  may  employ  two  models  for  the  extraction  of  Sn 
from  the  solder  drop  as  it  reacts  with  the  Cu  substrate.2  In  the  first  approach,  we  could  suppose  that 
the  boundary  behaves  like  a diffusion  couple.  In  such  a model,  the  concentration  at  the  substrate 
is  fixed  for  all  times  greater  than  the  initial  time, 


'l  • 


(2.6) 


For  example,  cL  may  be  taken  as  the  liquidus  concentration  from  the  equilibrium  phase  diagram  at 
the  temperature  of  the  substrate.2  A more  general  approach  would  be  to  impose  a mixed  transport 
condition  at  the  substrate, 

Dczi  + Hc(cl  — c ) = 0,  (2-7) 

where  the  coefficient  Hc  is  used  to  model  the  reaction  rate  at  the  substrate.  We  have  constructed  a 
model  using  boundary  condition  (2.7).  Equation  (2.6)  can  then  be  viewed  as  a limiting  case  of  (2. 
7)  for  large  Hc.  The  boundary  condition  (2.7)  is  a simple  model  for  the  dissolution  of  the  substrate; 
more  complicated  models  for  the  reaction  there  are  a subject  of  current  research.  We  assume  that 
the  substrate  remains  flat  and  immobile. 

On  the  free  surface  z'  = h'(x',t'),  we  have  that 


dti 


w — 


dt* 
n • Vc# 


dh’  , 

u , 

dx<  ’ 

0, 
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(2.8) 

(2.9) 


n • r'  ■ n 
t • T • n 


An  • V'T'  + y^(T‘  - Too) 

Oair 


K'a, 
t • w, 


0. 


(2.10) 

(2.11) 

(2.12) 


Here  T1  is  the  stress  tensor,  n and  t are  the  unit  normal  and  tangent  vectors  respectively,  with  n 
pointing  out  of  the  liquid.  The  heat  transfer  to  passive  gas  above  the  drop  is  modeled  by  the  ratio 
of  the  thermal  conductivity  of  the  gas  Aair  and  the  thermal  boundary  layer  thickness  6air.  We  also 
have  that  the  curvature 


K' 


(2.13) 


The  surface  tension  cr  is  given  by 


ct(T',c')  = aw-  7 t{T'  - Tw)  - 7 C{c'  - cj; 


(2.14) 


for  solder  with  the  given  form  for  <7,  7^  > 0 and  7C  < 0 when  considering  Sn  as  the  impurity  in  a 
Pb-Sn  drop. 

The  initial  conditions  are 


h\x\  0)  = h'0(x'),  c\x',  h0,  0)  = Cq(x'),  h'0[a0)  - 0,  a'(0)  = a0. 
At  the  contact  line,  initially  we  have 


dh[ 


dx 


y(a0)  = - tan0o 


(2.15) 


(2.16) 


where  60  is  the  initial  contact  angle  between  the  free  surface  of  the  drop  and  the  substrate.  At 
x'  = 0,  we  have  the  symmetry  conditions 


a1(0)  = £1(0)  = g(0)  = 0. 


dx 


dx' 


The  initial  volume  is  given  by 


ro.  0 

Vo  = 2 / tiQ(x')dx', 

Jo 

and  will  be  conserved  for  all  subsequent  times.  At  subsequent  times, 


a£,(-VV) 


tan  9(t '), 


ra  ^ ) 

Vo  = 2 / ti(x',t')dx\ 

Jo 


dh!  /N 

^7(°.  0 


d3h'  . dc'  . . 

(0,t')  = —(0,h’,t')  = 0. 


dx'  dx,3^~’~  ‘ dx1 

The  contact  line  is  assumed  to  move  according  to  the  modified  “Tanner  Law” 

a[,  = K[0(t>)  - eA]m- 


(2.17) 

(2.18) 

(2.19) 

(2.20) 
(2.21) 

(2.22) 
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here  k > 0 is  a constant  and  6 a is  the  advancing  equilibrium  contact  angle  for  a static  contact  line. 
Because  we  are  considering  only  drops  that  spread,  we  need  not  consider  contact  angle  hysteresis. 
As  pointed  out  in  Ref.  2,  the  constant  ac  in  the  modified  Tanner’s  Law  has  been  estimated  to  be 


a 

ac  = k0  — 


(2.23) 


where  ac0  = (7r/4)3/10  « 0.05  by  de  Gennes15  for  m = 3;  see  also  Ambrose  et  al.A  for  an  estimate  of 
k0  w 0.02,  based  also  on  the  work  of  de  Gennes.  The  magnitude  of  ac  stems  from  the  assumption 
that  the  dissipation  at  the  contact  line  results  from  the  competition  of  surface  tension  and  viscous 
forces.  A value  of  ac0  « 0.05  has  been  used  in  the  present  work  based  on  the  estimate  given  in 
Ref.  2;  hence  we  have  ac  = 1.35  X 103cm/s.  In  general,  however,  this  constant  should  be  determined 
for  each  material  system;28  for  the  case  of  silicone  oil  on  glass  the  mobility  ac  is  on  the  order  of 
10~3cm/s.  We  note  that  the  estimate  for  ac0  is  based  on  a perfectly  wetting  system;  we  have  made 
the  assumption  that  we  can  use  this  estimate  as  an  upper  bound  in  the  case  of  partial  wetting.  We 
have  also  made  the  assumption  that  6a  is  a constant  throughout  the  spreading  process;  as  pointed 
out  in  Ref  2,  it  need  not  be  true  that  6a  be  independent  of  all  other  variables  in  the  problem.  For 
example,  6a  may  decrease  during  spreading  because  the  action  of  the  solder  flux  may  modify  the 
substrate’s  surface  tension  during  spreading.  Determining  such  variation  and  incorporating  it  into 
a model  of  spreading  is  beyond  the  scope  of  this  paper.  Because  of  these  issues,  we  choose  to  deal 
with  only  the  planar  geometry  and  will  seek  qualitative  comparison  with  the  limited  spreading  data 
available  to  us. 

The  exponent  m for  the  modified  Tanner  law  is  also  a subject  of  ongoing  research;  a fit  of 
contact  angle  hysteresis  data  suggests27  m — 3,  while  recent  research  by  de  Gennes  and  coworkers 
suggests  a value  of  m — l.25  Experimental  determinations  have  found  m = 2.8  for  silicone  oil  on 
glass28  and  m = 2.6  for  silicone  oil  on  a silicon  wafer.37  We  shall  use  m — 3 in  this  work. 


3.  Lubrication  Theory  for  Reactive  Spreading 

We  next  derive  appropriate  governing  equations  for  the  drop  taking  into  account  the  transport  of 
both  solute  and  heat.  The  purpose  here  is  to  examine  the  effect  of  the  reaction  boundary  condition 
(2.7)  on  the  distribution  of  solute  and  the  subsequent  effect  on  the  spreading  of  the  drop.  This 
model  is  an  extension  of  the  one  considered  by  Ehrhard  and  Davis.27  The  boundary  conditions  for 
the  solute  are  different  from  the  temperature  field;  the  free  surface  allows  no  flux  of  solute,  while 
there  is  a mixed  boundary  condition  for  the  concentration  on  the  substrate.  For  the  thermal  field, 
we  retain  the  (mixed)  heat  transfer  condition  on  the  free  surface  and  a fixed  temperature  on  the 
substrate  as  in  Ref.  27. 

We  use  the  following  scales:27 


x = 


u = 


z — 


CLq 


u 


GO  00  5 


AC0J1  , 
t = — 
a0 


w 


w = 


ao0  o / 

V=^V 


T'  - T 

rp  CXI 

AT 


j 


c = ^,e  = fi,  v 


Ac 


0o 


W 

Vo 


(3.1) 
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Here  AT  = Tw  — T^,  Ac  = cm  — cL  and  cm  is  the  maximum  initial  concentration  in  the  drop. 
When  these  scalings  are  introduced  into  the  governing  equations,  we  obtain  at  leading  order  in  the 
limit  $o  — * 0, 


ux  + wz  = 0 (3.2) 

~Px  + Uzz  = 0,  (3.3) 

-?*-§  = 0>  (3-4) 

Tzz  = 0 (3.5) 

9qSR  (ct  + ucx  + wcz)  = czz  (3.6) 

where  the  subscripts  denote  partial  differentiation,  and  the  dimensionless  parameters  G,  C,R  and 
S will  be  discussed  shortly.  For  liquid  metals,  the  Schmidt  number  S = /x/ pD  is  typically  large; 
for  example,  5 = 81  for  Pb-Sn.  We  are  thus  motivated  to  retain  the  convective  terms  in  the  solute 
equation  in  anticipation  of  the  large  S balancing  0q  and  the  Reynolds  number  R = pht9™ao/ fi.  On 
2 — 0,  the  boundary  conditions  become 

(3.7) 

(3.8) 

(3.9) 


u 


Pc 


■u. 


- Hcc  = 0, 


w = 0. 


Here  we  have  taken  the  nondimensional  slip  coefficient  (3(h)  = (3a/ha]  we  shall  consider  only  a = 1 
in  this  work.  To  our  knowledge,  this  form  for  the  slip  was  first  used  by  Greenspan.18  On  2 = h(x,t), 
we  have 


w — ht  — uhx , 

(3.10) 

Cp  — hXX) 

(3.11) 

uz  — (cx  T hxcz)  — ^(Tx  + hxTz ), 

(3.12) 

c2  = 0, 

(3.13) 

Tz  + BT  = 0. 

(3.14) 

Initially,  we 

have  to  leading  order  that 

h(x,  0)  = h0(x),  h0(l)  = 0,  a(0)  = 1, 

(3.15) 

&»--§• 

(3.16) 

5>=  5>)  = o, 

(3.17) 

1 = 2 ho(x)dx, 

Jo 

(3.18) 

c = c(x,  0); 

(3.19) 

for  t > 0, 

h(a,  t)  = 0, 

(3.20) 
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dh 

dx 


(M)  = -0(t), 


dh  . \ d3h.  . dh 

= = = 0> 

,Q(0 


raV) 

1 = 2 h(x,t)dx. 

Jo 

The  contact  line  radius  is  given  by 

at{t ) = [0(2)  - 0^]m  . 


(3.21) 

(3.22) 

(3.23) 


(3.24) 


The  nondimensional  numbers  making  an  appearance  are 


C = — — Jr-  (Capillary  no.), 
G = (Bond  no.), 

0~W 


AC  = 


ACc 
Pi  = 


find™ 


7t(AT)0o 

fihz6™ 

7c(Ac)6>0 


(Thermocapillary  no.), 
(Solutocapillary  no.), 


Q' 

2l 2 (Shp  Coefficient). 
a090 


R = — -1 - (Reynolds  no.), 
NP 

S = (Schmidt  no.), 


B = 

nc  = 


Q’O^O  ^air  / dan 


A 

Hcao0o 

D 


(Biot  no.), 


(3.25) 

(3.26) 

(3.27) 

(3.28) 

(3.29) 

(3.30) 

(3.31) 

(3.32) 

(3.33) 


We  expect  that  the  product  RS  is  large  and  we  assume  that  the  Peclet  number  V = OqRS  is  0(1). 
The  nonlinear  convective  term  is  then  retained  in  the  solute  transport  equation. 

The  approximate  solution  will  follow  a mix  of  the  approaches  of  Ehrhard  and  Davis27  and 
Reisfeld  et  a/.;34  in  Ref.  34,  a Karman-Polhausen-like  approach  is  used  for  the  concentration  field 
in  a large  Peclet  number  limit  that  is  analogous  with  our  case.  The  idea  in  our  work  is  that 
the  profile  strongly  resembles  the  alloying  element  profiles  for  moderate  Peclet  number  transport  in 
drops.  This  assertion  is  based  on  finite  element  modeling  of  a drop  with  a fixed  contact  fine  position 
where  we  have  examined  the  development  of  the  flow  and  concentration  profiles  for  various  Schmidt 
numbers.36  The  linear  temperature  field  is  determined  independent  of  the  flow,  but  modifies  the 
flow  via  the  tangential  stress  condition. 

We  may  integrate  the  solute  equation  in  the  bulk  in  the  vertical  direction,  and  using  continuity, 
the  kinematic  condition  on  the  drop  surface,  and  the  flux  condition  on  the  substrate,  we  obtain 


d_ 

dt 


, d 

c dz  H 

I r\ 

OX 


[h 

/ uc  dz  = 
Jo 


Qz  (x , 0,  t) 

V~ 


(3.34) 
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(3.35) 


We  may  integrate  the  continuity  equation  in  a similar  manner  to  obtain 


ht  + 


dx  Jo 


udz  = 0. 


Ehrhard  and  Davis  were  able  to  solve  the  necessary  equations  in  the  non-reactive  case  and  derive 
a single  evolution  equation  for  the  free  surface  shape  via  equation  (3.35).  In  our  work,  the  mass 
transport  equation  (3.34)  is  not  easily  solved  and  so  we  are  motivated  to  find  an  approximate 
solution  for  the  solute  field  in  order  to  reduce  the  problem  to  evolution  equations  for  the  free 
surface  shape  and  the  appropriate  concentration  variable.  To  this  point  the  evolution  equations  (3. 
34)  and  (3.35)  are  independent  of  the  model  used  for  the  reaction  at  the  substrate. 

In  a manner  similar  to  Ref.  34,  we  pose  the  following  behavior.for  the  solute  field: 


c = c2(x,  t)  + Hcc2(x,t)z  + ca£(x,  t)zac, 


(3.36) 


where  ac  > 1 is  a constant;  this  form  is  consistent  with  the  flux  condition  (3.8)  on  the  substrate. 
In  order  to  satisfy  c2(x,/i,t)  = 0,  we  eliminate  cac  to  obtain 


c — c2  + 7tcc  2 


etc 


(3.37) 


Here  c2(x,t)  is  the  deviation  of  the  concentration  from  the  liquidus  on  the  substrate,  and  the 
flux  on  the  substrate  is  proportional  to  c2.  This  profile  resembles  those  seen  in  full  Navier- Stokes 
computations  for  a drop  with  a stationary  contact  line,  a moderate  contact  angle,  and  a moderate 
Peclet  number;36  this  profile  may  be  a poor  approximation  for  the  the  extremely  large  Peclet 
numbers  which  occur  in  an  actual  solder  drop.  It  corresponds  more  closely  to  a material  with  larger 
diffusivity  than  those  of  typical  liquid  metals.  Allowing  more  complicated  spatial  dependence  in 
the  Karman-Polhausen  profile  is  possible  at  the  expense  of  adding  more  free  parameters;  we  choose 
not  to  do  this. 

The  z-component  of  the  momentum  equation  (3.4)  can  again  be  used  to  find  the  pressure;  after 
using  the  normal  stress  boundary  condition  (3.11),  one  finds  that 

p(z,z,t)  = y(/i  - z)  - i/ixx.  (3.38) 


The  x-component  of  the  momentum  equation  (3.3)  can  then  be  used,  along  with  the  boundary 
conditions  to  find  that 


u{x,z,t)=  -i(D3x/l)y  + A2  (^2+  , 


where 


and 


A2  — — \ hD2xh  + M 


(1  + Bh)2 


- Mr 


C2x  + 1 ) He  ( c2h ) 


OLr 


D2xh  — (hxx  Gh)x. 

We  have  defined  the  thermal  and  solutal  Marangoni  numbers 

BC  C 

M = — , and  Mc  = — 

AC  ’ A Cc 


(3.39) 

(3.40) 

(3.41) 

(3.42) 
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We  may  integrate  u once  more  in  the  vertical  direction;  we  find  that 

r h ]_  f \ 

U = Jv  u dz  = ~-{D3xh)  — + A2\  — + Pij  • (3.43) 

Substitution  into  Eq.  (3.35),  ht  + ux  = 0,  gives  an  evolution  equation  for  the  shape  of  the  free 
surface  of  the  drop: 


ht  + — s ( — + (3\h  j D3xh-\- 


MK 


(i  + Bhy 


- Mr 


c2x  + (1  — “=■  ) hic(c2h)t 


= 0. 


(3.44) 


We  may  then  evaluate  the  integrals  required  in  Eq.  (3.34);  the  result  is  a second  evolution  equation, 
viz. 


where 


C2/1  + 7icbiC2h2 


-f  [C2U  -f  'HcC2A^]x 


(3.45) 


A4  — /1  + ^.2/25 


and 


h 


1 

2 c 


b3h4D3xh , 


f2  = b2h 3 + b1/31h, 


bk 


ac{ac  + k)  — (k  -f  1) 
ctc{ac  + k)(k  4-  1) 


k = 1,2,3. 


(3.46) 

(3.47) 

(3.48) 


The  evolution  equations  (3.44)  and  (3.45)  with  attendant  boundary  conditions  constitute  our  model 
for  reactive  spreading;  the  equations  in  the  limit  of  infinitely  large  7ic  may  be  considered  a second 
model  in  which  the  Dirichlet  condition  on  the  solute  field  is  c = 0.  Such  a model  would  consider 
the  reaction  at  the  substrate  to  be  classical  interdiffusion  with  infinitely  fast  kinetics.  We  have  also 
carried  out  computations  on  this  type  of  model;  those  results  are  not  reported  here. 


3.1  Boundary  and  initial  conditions 

The  evolution  equations  are  also  coupled  to  the  contact  line  motion  via 

at  = [~hx(a,t)  - Oa]™.  (3.49) 

Both  h and  c2  are  subject  to  the  symmetry  conditions  at  the  origin,  namely 

hx  — hxxx  — e2x  — 0*  (3.50) 

The  fluid  volume  is  also  conserved, 

1 /*a(<) 

-=  h(x,t)dx.  (3.51) 

2 Jo 
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This  is  an  approximation  in  our  model,  since  there  is  some  solute  leaving  the  drop,  but  we  assume 
that  the  volume  of  the  drop  is  not  affected  by  this.  For  the  initial  condition,  we  have  used 

c2(x,0)  = ehl(x),  (3.52) 

where  the  amplitude  e is  chosen  such  that  c(0,  h,  0)  = 1.  The  above  expression  for  the  initial  solute 
profile  must  satisfy  C2X{x  = a)  = 0]  this  boundary  condition,  in  addition  to  C2  = 0 are  seen  to  be 
natural  for  this  kind  of  problem  after  consideration  of  the  appropriate  wedge  flow  problems;40,41  in 
essence  we  have  matched  the  outer  solution  we  are  finding  onto  the  inner  wedge  flow  problem.  We 
thus  require  the  outer  solution  we  find  in  our  evolution  equation  to  have  the  correct  behavior  near 
the  contact  line. 

The  initial  drop  shape  for  all  results  presented  in  this  work  is  given  by 

h0(x)  = 5(1  - x2).  (3.53) 

We  have  adopted  this  initial  condition  for  the  purposes  of  comparison;  however,  because  we  have 
retained  the  time  derivative  in  the  evolution  equations,  we  could  have  used  any  reasonable  drop 
shape.  In  the  limit  of  small  capillary  number,  the  initial  distortions  from  a circular  cap  will  be 
eliminated  relatively  quickly  compared  to  the  time  scale  of  the  spreading. 

At  this  point  we  would  like  to  emphasize  that  these  initial  conditions  are  in  a sense  the  starting 
point  for  a long-time  solution  to  the  spreading  problem.  We  are  applying  this  kind  of  theory  for 
reactive  spreading  with  the  thought  that  the  initial  transients  involving  melting  are  completed  and 
fluid  dynamic  spreading  of  a roughly  circular  cap  is  beginning. 


4.  Numerical  Integration  of  the  Evolution  Equations 


The  above  equations 
then  have 

and 


were  mapped  onto  a fixed  domain  — 1 < f < 1 by  defining  ( = x/a{t );  we 
dh(x,t)  1 dh(£,t) 


dx 


dh(x,t) 

dt 


d_ 

dt 


2 

Off  _d_ 
a d(. 


(4.1) 

(4.2) 


The  resulting  equations  are  then  discretized  pseudospectrally  using  Chebyshev  polynomials  and 
collocated  at  the  Gauss-Lobatto  points.38  An  advantage  of  this  approach  is  that  we  are  able  to  work 
in  the  physical  domain  and  use  a differentiation  matrix  to  evaluate  the  derivatives  at  the  collocation 
points  while  obtaining  high  order  accuracy  in  the  discretization.  This  results  in  a system  of  ODEs 
in  time  for  the  elevations  of  the  drop  surface  at  the  collocation  points  and  the  radius  of  the  drop, 
as  well  as  the  algebraic  condition  which  conserves  volume.  Due  to  the  symmetry  of  the  drop,  we 
need  only  solve  the  equations  on  0<f  <i.  This  system  of  differential- algebraic  equations  (DAEs) 
is  then  solved  by  the  package  DASSL,39  which  uses  an  adaptive,  implicit  time  stepping  scheme 
that  is  up  to  5th-order  accurate.  Our  solution  to  the  DAE’s  is  a more  accurate  approach  than  the 
first  order  scheme  (in  time)  used  by  Haley  and  Miksis;35  otherwise  the  approaches  are  practically 
identical. 
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5.  Results 


The  parameters  for  Pb-Sn  solder  near  the  eutectic  point  as  well  as  initial  data  used  in  the  calcu- 
lations are  given  in  Table  1.  We  point  out  that  these  parameters  are  what  we  would  like  to  use 
for  computation  for  solder  drops.  We  are  able  to  achieve  the  desired  parameters  for  a non-reactive 
drop  in  Section  5.1;  for  the  reactive  case  we  discuss  what  we  can  achieve  in  Section  5.2. 

5.1  Non-reactive  Spreading 

For  this  case,  the  drop  is  at  uniform  concentration.  The  evolution  equations  reduce  to  a single 
evolution  equation  for  the  free  surface  shape;  that  equation  is  a modified  version  of  Ehrhard  and 
Davis’  result,  viz., 

= 0.  (5.1) 

The  modification  to  their  model  is  that  we  solve  the  evolution  equation  while  retaining  the  time 
derivative  and  a nonzero  Biot  number  in  the  denominator  of  the  Marangoni  term.  This  equation  is 
subject  to  the  symmetry  conditions  at  the  origin,  the  contact  line  constitutive  law,  and  conservation 
of  volume. 

In  Figure  1,  we  show  results  for  the  pure  capillary  case  G = M = 0;  numerical  results  using 
the  technique  described  above  are  shown  for  two  different  advancing  equilibrium  contact  angles, 
0,4  = 0.1  and  0,4  = 0;  Figure  2 shows  the  (approximate)  slope  S of  the  In  a vs.  In  2 plot.  Since  for 
8 =constant  the  growth  of  the  radius  is  proportional  to  t6 , we  refer  to  6 as  the  spreading  exponent. 
The  chain-dashed  curve  shown  in  these  figures  is  the  solution  to  the  C — 0 problem.15, 18,27  In  the 
planar  geometry,  for  C — 0 we  must  solve  an  ODE  for  the  contact  line  radius 

at  = (2?  - ®-4)  ’ a(0)  = L (5-2) 

The  steady-state  contact  line  radius  is  given  by  a ^ — ^3/ (20^).  Good  agreement  at  long  times  is 
obtained  with  the  numerical  method;  it  is  expected  that  the  two  should  agree  at  long  times  since 
the  freedom  to  specify  an  arbitrary  initial  shape  for  the  drop  shape  is  given  up  when  C = 0 in  the 
problem.  The  final  radius  is  attained  very  accurately.  We  also  note  that  the  radius  of  the  drop 
increases  in  proportion  to  t6 , 8 = 1/7  for  0^  = 0 in  Figure  2,  as  is  expected  for  the  planar  geometry 
and  capillary  spreading.24,27  We  view  this  case  as  a verification  of  the  numerical  method,  especially 
for  the  long  times  required  in  the  computation.  It  is  clear  that  t 1//7  spreading  is  not  obtained  for 
nonzero  advancing  equilibrium  contact  angles. 

We  also  note  that  in  some  work,4,5  the  C = 0,4  = 0 (perfect  wetting)  solution  has  been  applied 
to  the  spreading  of  solder  droplets  where  0^  ^ 0 (partial  wetting).  It  is  clear  why  the  discrepancy 
in  time  scale  between  the  experiments  and  the  theory  exists  in  light  of  the  difference  exhibited  in 
Figures  1 and  2. 

The  effect  of  gravity  is  included  in  results  shown  in  Figures  3 and  4.  In  these  figures,  cases  for 
0,4  = 0 and  0^  = 0.1  are  shown  as  well  as  a quasistatic  case  C — 0,  0,4  = 0.1.  The  simplified 
problem  that  is  solved  for  C = 0 is  given  by  (after  Ref.  27) 

G 

a‘  ~ ["2(1  - aG1/2cothaG1/2)  " 


, a(0)  = 1.  (5.3) 


. 1 
ht  + c 


h3  a 

Y + Plh 
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Long  time  agreement  is  again  observed  between  the  numerics  and  the  quasistatic  solutions  of  Ref.  27. 
The  C — 0 solution  spreads  faster  at  early  times,  and  then  the  two  become  indistinguishable  at 
late  times;  see  also  Table  2.  At  long  times,  we  recover  t1/4  spreading  in  the  0^  = 0 case;  this 
gravity-dominated  regime  occurs  after  an  early  transient  which  resembles  pure  capillary  spreading 
(t1/7).  This  behavior  is  what  is  expected  based  on  previous  experimental  and  theoretical  work.24,27 
Once  again,  power-law  spreading  is  not  achieved  for  0^  ^ 0. 

Thermal  Marangoni  effects  are  included  in  the  results  shown  in  Figures  5 and  6;  here  we  have 
neglected  the  effect  of  gravity.  In  these  figures,  the  thermal  Marangoni  effect  retards  spreading, 
as  would  be  expected  in  solder  drops  that  are  heated  from  below.  We  find  that  the  final  radius 
agrees  well  with  the  asymptotic  result  given  in  Ref.  27  for  M C 1 as  given  in  Table  3.  The  data 
for  Table  3 is  for  0^  = 0.5.  The  magnitude  of  the  Marangoni  effect  given  by  M — 0.01  is  expected 
to  be  the  maximum  which  could  be  reasonably  expected  in  a solder  drop  under  the  conditions  of 
the  experiments  in  Refs.  1,2  and  is  based  on  the  data  in  Refs.  6,7.  For  0^  = 0.1,  the  computed 
final  radius  is  reduced  by  about  10%  from  the  M — 0 case;  this  is  a much  bigger  effect  than  for  the 
larger  0^  data  in  the  table. 

As  is  shown  in  Figure  6,  there  is  no  longer  true  power-law  spreading  for  0^  ^ 0,  but  there  is  a 
period  at  early  times,  where  a roughly  constant  6 exists.  A rapid  decrease  in  the  rate  of  spreading 
occurs  at  later  times,  and  our  numerical  solution  suggests  that  spreading  stops  even  for  0^  = 0 
as  predicted  by  Ref.  27.  Our  numerics  break  down  before  the  final  radius  is  reached  however,  due 
to  the  development  of  singularities  in  the  higher  derivatives,  as  pointed  out  in  Ref.  27.  In  order 
to  accurately  represent  the  approach  to  the  final  radius  without  enforcing  a quasistatic  behavior 
would  require  more  sophisticated  numerical  methods  than  we  have  employed  here. 

Figure  7 shows  radius  vs.  time  results  when  gravity  and  thermocapillarity  are  included.  When 
the  Marangoni  number  is  positive  (negative),  the  temperature  field  corresponds  to  a hotter  (colder) 
substrate  than  the  surroundings;  this  situation  retards  (promotes)  spreading  compared  to  the  M = 
0 case  when  gravity  is  included.  For  the  parameters  used  here,  the  case  when  Marangoni  flow 
retards  spreading  with  G = 0.15  is  very  close  to  the  case  M = G = 0 shown  in  Figure  1. 

In  Figure  8,  we  display  results  that  show  the  sensitivity  of  the  model  to  the  slip  parameter  (3\ 
when  the  Marangoni  effect  is  included.  Greenspan18  showed  that  when  the  capillary  number  is 
small,  the  solutions  for  the  spreading  drop  are  insensitive  to  the  slip  parameter;  in  the  absence  of 
Marangoni  effects,  our  numerics  bear  this  out.  When  Marangoni  effects  are  included  however,  we 
must  have  a small  slip  parameter  in  order  to  achieve  the  correct  final  radius. 

In  Figure  9,  we  show  the  shape  of  the  free  surface  of  the  drop  including  gravity  and  thermo- 
capillarity for  an  advancing  contact  angle  0^  = 0.  The  drops  are  thin  in  physical  coordinates,  but 
do  not  necessarily  appear  thin  in  the  scaled  coordinates.  At  late  times,  an  inflection  point  appears 
near  the  outside  edge  of  the  drop;  in  this  case,  the  radius  at  the  end  of  the  computation  is  almost 
7.5.  The  case  shown  for  0,4  = 0 displays  this  inflection  most  clearly.  For  nonzero  0^,  the  inflection 
point  is  difficult  to  detect  visually. 

5.2  Reactive  Spreading 

Figures  10  and  11  show  results  from  the  numerical  solution  of  the  evolution  equations  (3.44)  and 
(3.45);  here,  V = 200,  7tc  — 2000,  C - .05,  0^  = 0.1,  G — 0 and  M — —Mc  = 0.01  which  we 
designate  case  I.  We  may  think  of  the  spreading  as  taking  place  between  the  angles  of  45°  and  3°; 
then  t = 6600  is  roughly  one  second.  This  choice  of  0^  implies  a 15  times  reduction  in  the  slope  of 
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the  drop  surface  at  the  contact  line  during  spreading.  We  have  computed  solutions  in  this  section 
for  drops  whose  volume  has  initial  radius  a0  = 0.17cm  at  90  = 45°.  The  Marangoni  numbers  were 
hxed  at  there  values  based  on  the  estimates  in  Table  1,  with  the  expectation  that  the  two  effects 
are  of  comparable  magnitude. 

We  see  in  Figure  10  that  at  early  times  the  solutal  Marangoni  effect  retards  spreading  since 
the  radius  at  each  instant  of  time  lags  that  of  the  corresponding  conditions  for  the  non-reactive 
drop  with  or  without  the  thermal  Marangoni  effect.  This  slowing  of  the  spreading  relative  to  non- 
reactive cases  disappears  relatively  quickly;  the  excess  solute  in  the  drop  has  been  extracted  by  the 
substrate  and  the  slowing  of  the  contact  line  abates.  The  absence  of  solute  (i.e.,  the  arrival  of  c!  at 
cL ) then  eliminates  the  solutal  Marangoni  effect. 

In  Figure  11,  the  concentration  of  solute  at  the  top  of  the  drop  and  the  normalized  total  amount 
of  solute  are  shown  as  a function  of  time.  It  is  clear  that  the  drop  has  an  increased  concentration 
at  the  top  of  the  drop  when  compared  to  the  initial  state  for  some  times.  The  total  amount  of 
solute  in  half  of  the  drop  Ic{t ) decays  monotonically,  as  shown  by  the  lower  curve,  where 

ra  rh 

Ic(t)  = / / c(x,z,t)dxdz  (5-4) 

Jo  Jo 

The  increase  in  the  concentration  at  the  top  of  the  drop  is  an  artifact  of  the  model,  since  diffu- 
sion would  not  allow  such  an  increase;  however,  since  the  increase  is  not  large,  we  are  willing  to 
accept  such  behavior  under  our  stated  philosophy  of  giving  Marangoni  effects  the  most  favorable 
circumstances  to  impact  spreading.  It  is  perhaps  not  surprising  that  such  behavior  is  possible  in 
our  model  because  we  are  not  solving  the  mass  transport  equation  pointwise,  but  we  are  really 
satisfying  it  in  an  integrated  sense  by  using  a Karman-Polhausen  profile  for  the  solute  field.  The 
increase  in  local  concentration  does  not  appear  in  problems  where  there  is  a fixed  concentration  on 
the  moving  boundary  as  in  Reisfeld  et  al.  33,34 

For  completeness  we  show  the  spreading  exponent  behavior  in  Figure  12.  Initially,  the  spreading 
is  slowed  and  6 is  reduced,  but  it  increases  later  in  the  spreading  while  the  spreading  catches  up  to 
the  thermal  case.  At  early  times,  spreading  is  retarded  by  the  presence  of  the  solute  field,  and  late 
in  the  spreading  the  thermal  field  retards  the  spreading. 

The  behavior  of  the  solute  field  and  stream  function  are  shown  in  detail  in  the  sequence  depicted 
in  Figures  13  through  17.  The  top  (bottom)  part  of  the  figure  shows  the  concentration  field  (stream 
function).  By  a fairly  early  time,  the  solute  concentration  has  become  very  small  (below  the  error 
tolerance  of  the  time  integration)  and  by  that  time  the  solutal  Marangoni  effect  plays  no  further 
role.  In  the  time  that  it  is  operating  however,  it  generates  a vortex  which  retards  the  spreading  in 
comparison  with  the  case  with  no  Marangoni  effects  (cf.  Figure  12).  It  is  important  to  note  that 
the  maximum  of  the  stream  function  is  near  10_1  initially  and  is  about  10~5  at  the  latest  time  of 
the  calculation.  This  wide  variation  in  the  strength  of  the  flow  and  the  fact  that  the  third  derivative 
of  the  free  surface  shape  is  required  to  compute  the  stream  function  makes  this  computation  an 
extreme  test  for  the  collocation  method  we  have  used.  Some  error  due  to  aliasing  is  visible  in  the 
streamlines  at  the  latest  times;  the  numerical  error  is  compounded  by  interpolation  of  the  solution 
from  the  Chebyshev  points  (32  or  48  on  the  half  interval)  to  a fine  uniform  mesh  for  contour  plotting 
(200  points).  The  time  scale  for  the  exit  of  the  solute  from  the  drop  is  set  by  the  Peclet  number 
V and  to  a lesser  extent  by  7ic\  the  time  scale  of  the  spreading  flow  is  set  primarily  by  k,  m , and 
0^.  We  have  used  the  largest  V for  which  we  could  compute  solutions.  Based  on  the  diffusivity  in 
Table  1,  a Peclet  number  of  V = 104  to  106  would  occur  in  the  actual  drop.  For  values  of  V < 100, 
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a vortex  due  to  the  solutal  Marangoni  effect  did  not  form. 

In  carrying  out  the  lubrication  approximation,  an  approximate  no-flux  boundary  condition  was 
applied  at  the  free  surface,  namely  cz(x,h,t)  = 0.  Mass  flux  through  the  free  surface  of  the  drop, 
given  by 


/ 


j h , t)dx, 


(5.5) 


compared  to  that  through  the  base, 

f cz(x,0,t)dx , (5.6) 

Jo 

is  at  most  about  0.3#o  at  t = 0,  and  rapidly  decreases  to  less  than  0.05#o-  Clearly  the  approximation 
of  the  no  flux  boundary  condition  improves  as  the  initial  contact  angle  do  decreases. 

An  example  of  spreading  with  less  extreme  deformation  of  the  drop  which  includes  gravity  and 
Marangoni  effects  is  shown  in  Figures  18  through  22,  which  we  call  Case  II.  In  this  case,  V = 200, 
7ic  — 2000,  C — 0.05,  0,4  = 0.5,  G = 0.15  and  M — —Mc  — 0.01.  Here  we  could  think  of  the 
spreading  as  taking  place  between  the  angles  of  30  and  10  degrees;  then  t = 1900  is  roughly  one 
second.  The  solutal  Marangoni  effect  dies  off  during  the  time  the  thermal  Marangoni  effect  becomes 
important.  If  less  spreading  occurs,  the  time  scales  of  the  solutal  and  thermal  Marangoni  effects 
overlap  for  a wider  range  of  parameters,  which  results  in  (in  this  case)  “addition”  of  the  two  effects. 
While  this  particular  run  illustrates  clearly  the  effects  of  fluid  dynamics,  it  is  relatively  brief  in 
terms  of  the  duration  of  spread  for  typical  solder  spreading. 


6.  Discussion  and  Conclusion 

The  computational  approach  to  the  spreading  of  a drop  was  verified  against  previously  known 
solutions  in  the  presence  of  capillarity,  gravity,  and  thermocapillarity  (Marangoni  effect)  given  by 
Ehrhard  and  Davis.27  The  good  agreement  of  our  results  for  the  non-reactive  cases  indicates  that 
the  method  works  reasonably  well  for  nonzero  advancing  equilibrium  contact  angles.  When  the 
advancing  equilibrium  contact  angle  is  zero  (perfect  wetting  case),  we  find  that  our  numerical 
approach  breaks  down  late  in  the  spreading  process  when  the  Marangoni  effect  is  present.  Because 
the  rate  of  change  of  the  radius  is  slow  by  this  time,  we  expect  that  the  computed  answer  is 
a reasonable  estimate  of  the  final  radius.  As  expected,  the  agreement  between  our  computed 
solutions  and  the  C = 0 solutions  improved  with  increasing  time.  The  agreement  improves  for 
decreasing  C as  well. 

An  interesting  feature  from  our  computational  results  is  the  dependence  of  the  final  radius 
on  the  slip  coefficient  when  the  Marangoni  effect  is  present.  For  the  parameters  we  used,  a slip 
coefficient  less  than  about  (3\  = 10-4  was  required  to  get  a final  radius  that  was  insensitive  to  the 
value  of  the  slip  coefficient. 

In  contrast  with  the  drops  of  uniform  concentration,  the  results  of  the  computations  in  the  reac- 
tive case  have  no  limiting  analytical  solution  for  comparison.  The  solute  field  generates  complicated 
flow  patterns  early  in  the  computation;  it  may  generate  a vortex  within  the  flow  field  that  would 
otherwise  be  dominated  by  capillarity  when  the  solute  held  is  absent.  Depending  on  the  extent 
of  spreading,  the  model  may  allow  for  an  intermediate  period  where  capillarity  dominates  before 
thermal  Marangoni  effects  take  over.  Thus  the  results  also  indicate  that  the  solute  diffuses  out  to 
the  substrate  relatively  quickly  compared  to  the  spreading  of  the  drop.  The  theory  presented  here 
with  reactive  effects  corresponds  to  a fluid  that  has  higher  diffusivity  than  Pb-Sn  solder,  and  this 
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is  certainly  a large  part  of  the  cause  of  the  relatively  fast  behavior  of  the  solute  field.  The  higher 
diffusivity  used  in  the  results  of  this  paper  is  apparent  upon  solving  for  the  diffusivity  given  the 
Peclet  number  V = 200.  It  is  then  clear  we  are  operating  in  the  range  from  about  10-3  to  about 
10_1cm2/s  for  the  diffusivity,  depending  on  the  initial  angle  and  radius  given  the  estimate  we  use 
for  k.  It  may  possible  to  build  more  complex  behavior  into  the  spatial  shape  of  the  solute  field,  and 
presumably  to  access  more  extreme  diffusivities,  but  this  comes  at  the  cost  of  adding  additional  free 
parameters  into  the  problem.  The  possibility  remains  that  the  our  model  may  apply  to  situations 
other  than  metal/metal  systems,  e.g.,  organics. 

We  have  proposed  a model  for  reactive  partial  wetting  situation  that  attempts  to  characterize 
spreading  in  many  metal/metal  systems.  One  of  the  central  features  in  the  model  we  propose  is 
to  incorporate  a constitutive  law  for  the  contact  line  motion  in  order  to  treat  a macroscopic  fluid 
dynamics  problem.27  Use  of  the  assumption  of  perfect  wetting  ( 9 a = 0)  with  this  constitutive  law 
leads  to  orders-of-magnitude  discrepancy  with  experimental  results  as  pointed  out  in  Yost  et  a/.5 
and  Ambrose  et  al.4  While  the  geometry  of  the  model  is  planar,  rather  than  axisymmetric,  we  find 
that  the  time  scale  to  complete  nearly  all  of  the  spreading  is  definitely  of  the  correct  order  (about 
1-10  seconds  depending  on  the  flux  used  for  some  unpublished  spreading  results  for  fluxed  Pb-Sn 
on  several  different  Cu  and  Cu-Sn-intermetallic  substrates42).  The  same  can  be  said  for  the  time 
scales  of  previous  theories  of  Boettinger  et  al .,2  (although  they  did  not  point  this  out  in  their  work) 
and  of  Ehrhard  and  Davis.27 

Many  issues  remain  to  be  addressed  in  this  area.  The  contact  angle  may  depend  on  many 
variables  and  parameters  in  the  problem.2  We  have  also  neglected  the  effects  of  flux,  which  are 
important  in  terms  of  wettability,  but  may  also  be  important  in  terms  of  the  fluid  dynamics  during 
wetting.  A more  complete  treatment  of  the  large  Peclet  number  transport  would  be  desirable. 
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Density 

Viscosity 

Diffusion  coefficient 
Surface  tension 
Thermal  Marangoni  effect 
Solutal  Marangoni  effect 
Thermal  conductivity 
Heat  capacity 
Temperature  difference 
Concentration  difference 
Contact  line  mobility 


p 8.3  gm/cm3 
p 1.8x  10-2  g / (cm  s) 

D 2.7x  10~5  cm2/s 
cfw  500  erg/cm2 
7 t 0.05  erg/(cm2  K) 

7C  -1.0  erg/(cm2  wt  %) 

A 5.0xl06  erg/(cm  K S) 
cp  2.2  xlO6  erg/(gm  k) 
AT  100  K 
Ac  1.0  wt  % 
k,  1.35  xlO3  cm/s 


Table  1:  The  material  parameters  are  for  Pb-Sn.  The  parameters  for  the  passive  gas  (air)  above 
are  assumed  to  be  such  that  the  Biot  number  is  fixed  at  B = 0.1.  The  contact  line  mobility  and 
the  temperature  and  concentration  differences  are  estimates. 


G 

Asymptotic 

Numerical 

0.001 

1.732224 

1.732224 

0.01 

1.733783 

1.733784 

0.1 

1.749371 

1.749481 

1.0 

1.905256 

1.915288 

Table  2:  A comparison  of  our  numerical  results  with  the  asymptotic  results  of  Ehrhard  and 
Davis27  for  the  final  radius  a ^ of  an  isothermal  drop  with  gravity.  Here  0^  = 1/2,  and  so 
aP0 0 = VW07)  = \/3.  Their  asymptotic  results  are  given  by  a ^ = u^ofl  + (a^o)26?/30]. 
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Radius  a(t ) for  pure  capillary  spreading.  The  C — 0 case  is  for  the  spreading  of  a circular 
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Figure  2:  Spreading  exponent  8 for  pure  capillary  spreading.  The  C — 0 case  is  for  the  spreading 
of  a circular  cap. 


22 


a(  t) 
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Figure  3:  Radius  a(t ) for  spreading  with  capillarity  and  gravity.  The  C — 0 case  is  for  quasistatic 
spreading  of  a cap  with  gravity. 
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Figure  4:  Spreading  exponent  8 for  M = 0.  The  C — 0 case  is  for  quasistatic  spreading  of  a cap 
with  gravity.  The  spreading  exponent  asymptotes  to  6 = 0.25  for  the  perfectly  wetting  case. 


24 


Figure  5:  Radius  a{t ) for  spreading  with  thermocapillarity. 
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Figure  6:  Spreading  exponent  6 for  capillary  spreading  with  thermocapillarity. 
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Figure  7:  Radius  a(t)  for  capillary  spreading  with  gravity  and  thermocapillarity.  For  M — 0.01, 
the  Marangoni  effect  retards  spreading;  the  case  with  the  opposite  sign  promotes  spreading. 
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a ( t ) 
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Figure  8:  Radius  a(t)  for  capillary  spreading  with  gravity  and  thermocapillarity  as  a function  of 
(3\.  The  small  capillary  number  solutions  now  depend  on  the  slip  coefficient. 
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Figure  9:  Droplet  shapes  for  various  times  for  capillary  spreading  with  gravity  and  thermocapillar- 
ity. 
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Figure  10:  Drop  radius  as  a function  of  time.  Only  thermo-  and  soluto- capillary  effects  are  included 
here.  At  longer  times,  the  two  curves  involving  the  Marangoni  effects  merge. 
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Figure  11:  Concentration  at  the  top  of  the  drop  as  a function  of  time,  and  the  total  amount  of  solute 
Ic(t)  normalized  by  its  initial  value  Jc(0).  Only  thermo-  and  soluto- capillary  effects  are  included 
here. 
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Figure  12:  Exponent  of  drop  radius  as  a function  of  time.  Only  thermo-  and  soluto-capillary  effects 
are  included  here.  The  slowing  of  the  growth  of  the  radius  at  early  times  can  be  seen  by  a lowered 
6,  but  at  later  times  the  radius  catches  up  to  the  case  of  the  pure  drop  with  M = 0.01. 
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M 

Asymptotic 

B = 0 

B = 0.1 

0.001 

1.73055 

1.73067 

1.72943 

0.01 

1.71705 

1.71859 

1.71907 

0.1 

1.58205 

1.61122 

1.61539 

0.25 

1.35705 

1.47430 

1.48254 

Table  3:  A comparison  of  our  numerical  results  with  the  asymptotic  results  of  Ehrhard  and  Davis27 
for  the  final  radius  of  an  nonisothermal  drop  without  gravity.  The  columns  B = 0 denotes  the 
case  when  M / 0,  but  B = 0 in  the  denominator  of  the  Marangoni  term  in  the  evolution  equation 
(for  direct  comparison  with  Ref.  27).  Here  0^  = 1/2,  and  = y/3.  Their  asymptotic  results  are 
given  by  a = a^[l  + Ma^/(40^)]. 
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Figure  13:  Case  I,  t = 0.  The  top  figure  is  the  initial  dimensionless  concentration  field  c (in 
steps  of  Ac  = 0.1)  and  the  bottom  is  the  initial  stream  function  ip  (in  steps  of  Aip  = 0.01).  The 
concentration  and  stream  function  are  both  zero  at  the  substrate  for  this  and  all  subsequent  figures. 
The  stream  function  indicates  the  dominance  of  capillary  spreading  initially. 
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Figure  14:  Case  I,  t = 3.0,  Ac  = 0.1,  and  Ai p = 0.001.  The  stream  function  is  already  being 
strongly  affected  by  the  solute  gradient  on  the  free  surface,  and  the  concentration  at  the  top,  center 
of  the  drop  has  increased  due  to  the  flow  sweeping  the  alloying  element  there. 
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Figure  15:  Case  I,  t = 22.0,  Ac  = 0.1,  and  At/?  = 0.0002.  A vortex  has  formed  as  a result  of 
the  solute  gradient.  Note  the  rapid  decrease  of  the  maxima  of  the  stream  function  and  the  solute 
profile. 
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Figure  16:  Case  I,  t = 62.9  and  Aip  = 0.0001.  The  stream  function  is  no  longer  affected  by  the 
concentration  since  the  maximum  concentration  is  now  well  below  1%  of  the  initial  value;  we  no 
longer  show  the  concentration  field  because  of  this.  Note  the  continued  decrease  in  the  maximum 
of  the  stream  function.  The  stream  lines  do  not  indicate  strong  effects  from  thermocapillarity  yet. 
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Figure  17:  Case  I,  t = 7071.,  and  Aip  = 0.000005.  The  maximum  concentration  field  has  decreased 
to  far  below  1%  of  its  original  maximum,  and  so  we  no  longer  display  its  contours.  A single  vortex 
which  is  rather  weak  in  magnitude  occupies  the  drop.  The  spreading  of  the  drop  is  about  86% 
complete. 
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Figure  18:  Case  II,  t — 0.  The  top  figure  is  the  initial  dimensionless  concentration  held  c (in  steps 
of  Ac  = 0.1)  and  the  bottom  is  the  initial  stream  function  t/’  (in  steps  of  A ip  = 0.01)  for  case  II. 
The  stream  function  indicates  the  dominance  of  capillary  spreading  initially;  the  streamlines  have 
been  turned  somewhat  by  the  presence  of  gravity. 
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Figure  19:  Case  II,  t = 3.0,  Ac  = 0.1,  and  A ip  = 0.001.  The  stream  function  is  already  being 
strongly  affected  by  the  solute  gradient  on  the  free  surface,  and  the  concentration  at  the  top  of  the 
drop  has  increased  due  to  the  flow  sweeping  the  alloying  element  there. 
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Figure  20:  Case  II, 
the  solute  gradient, 
profile. 


t = 22.0,  Ac  = 0.1,  and  Aip  = 0.0002.  A vortex  has  formed  as  a result  of 
Note  the  rapid  decrease  of  the  maxima  of  the  stream  function  and  the  solute 
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Figure  21:  Case  II,  t = 62.9,  Ac  = 0.01,  and  A jp  = 0.0001.  Note  the  smaller  step  size  between 
concentration  contours  and  the  continued  decrease  in  the  maximum  of  the  stream  function.  The 
stream  lines  do  not  indicate  strong  effects  from  thermocapillarity  yet. 
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Figure  22:  Case  II,  t = 1465.,  and  A ip  = 0.00005.  The  maximum  concentration  field  has  decreased 
to  far  below  1%  of  its  original  maximum,  and  so  we  no  longer  display  its  contours.  A single  weak 
vortex  occupies  the  drop.  The  spreading  of  the  drop  is  about  98%  complete. 
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